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Abstract 

The Iberian Peninsula is a region with a high endemicity of species of the terrestrial slug subgenus Mesa- 
rion. Many of these species have been described mainly on subtle differences in their proximal genitalia. It 
therefore remains to be investigated 1) whether these locally diverged taxa also represent different species 
under a phylogenetic species concept as has been shown for other Mesarion species outside the Iberian 
Peninsula, and 2) how these taxa are phylogenetically related. Here, we analysed DNA sequence data of 
two mitochondrial (COI and 16S) genes, and of the nuclear ITSl region, to explore the phylogenetic 
affinities of two of these endemic taxa, viz. Arion gilvus Torres Mfnguez, 1925 and A. ponsi Quintana 
Cardona, 2007. We also evaluated the use of these DNA sequence data as DNA barcodes for both species. 
Our results showed that ITS did not allow to differentiate among most of the Mesarion molecular opera- 
tional taxonomic units (MOTUs) / morphospecies in Mesarion. Yet, the overall mean p-distance among 
the Mesarion MOTUs / morphospecies for both mtDNA fragments (16.7% for COI, 13% for 16S) was 
comparable to that between^, ponsi and its closest relative ^4. molinae (COI: 14.2%; 16S: 16.2%) and to 
that between yl. ^i&Mj and its closest relative yl. urbiae {COI: 14.4%; 16S: 13.4%). Hence, with respect to 
mtDNA divergence, both A. ponsi and A. gilvus, behave as other Mesarion species or putative species-level 
MOTUs and thus are confirmed as distinct 'species'. 
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Introduction 

The genus Arion Ferussac, 1819 is the most species rich genus of the terrestrial slug 
family Arionidae (Mollusca, Pulmonata, Gastropoda). It comprizes approximately 40 
species, grouped into four subgenera, viz. Arion s.s. Ferussac, 1819, Kobeltia Seibert, 
1873, Carinarion Hesse, 1926 and Mesarion Hesse, 1926. Species of the subgenus 
Mesarion (type species: Limax subfuscus Draparnaud, 1805) are characterized by 1) a 
medium body-size (up to 75 mm when extended), 2) an orange to dark brown dor- 
sum, 3) two dark bands on the sides of the mantle, 4) (usually) yellow to orange body 
mucus, and 5) an enlarged free-oviduct with a long and V-shaped ligula (Kerney et 
al. 1983). M^nj Mesarion species are highly polymorphic with respect to body colour 
and genital anatomy. As a consequence, the species limits and phylogenetic relation- 
ships of taxa within this subgenus have been debated for decades (e.g. Garrido et al. 
1995, Castillejo 1997, 1998, Pinceel et al. 2004, 2005a, b, Quinteiro et al. 2005). 
Arion subfuscus (Draparnaud, 1805) (type locality: Montagne Noire, France) is prob- 
ably the most problematic "species" within Mesarion as it shows an overwhelming 
amount of variation in body pigmentation, genital anatomy, and reproductive behav- 
ior [see Garrido et al. (1995) and the references listed in their table 1]. This variation 
often has been interpreted as indicating reproductive isolation between geographi- 
cally isolated populations, and A subfuscus thus is considered a species complex (Wik- 
tor 1973, Walden 1976, De Winter 1986, Backeljau 1989, Altonaga et al. 1994, 
Backeljau et al. 1994, Garrido et al. 1995). Especially in the Pyrenees and the coastal 
regions of Spain there are local, morphologically diverged populations (e.g. Garrido 
et al. 1995, Castillejo 1998). Several of these have been described as endemic species 
on the basis of where the epiphallus, oviduct and pedunculus of the bursa copulatrix 
open into the atrium, in combination with differences in the relative lengths of the 
vas deferens and the epiphallus (e.g. Castillejo 1998, Garrido et al. 1995, Quintana 
Cardona 2007). Two of these endemic taxa occur in the eastern coastal region of 
Spain or the Balearic Islands, viz. Arion gilvus Torres Mi'nguez, 1925 and A. ponsi 
Quintana Cardona, 2007. 

Arion ponsi (Figure 1) was described from Menorca (Balearic Islands, type local- 
ity: Barranc dAlgendar). The species has a medium body size (range: 54-66 mm), an 
orange to beige dorsal body colour with dark lateral bands that can be blurry in the 
posterior parts, a foot sole that is cream coloured with a greyish hue, and a transparent 
body mucus (Quintana Cardona 2007). Its genital anatomy is very similar to that of^. 
gilvus, A. iratii Garrido, Castillejo & Iglesias, 1995, A molinae Garrido, Castillejo & 
Iglesias, 1995 and^. lizarrustii Garrido, Castillejo & Iglesias, 1995, but its epiphallus 
is shorter than the vas deferens (as in A. molinae) and opens into the genital atrium in 
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Figure I . Arion ponsi Quintana Cardona, 2007 from Menorca (Balearic Islands, Spain). 

between the oviduct and the pedunculus of the bursa copulatrix (unUke in A. molinae, 
where the pedunculus is positioned in between the epiphallus and oviduct) (figures 
3-5 in Quintana Cardona 2007). 

Arion gilvus (Figure 2) was described from 'Mandol' in the Spanish Province of 
Tarragona. However, the toponym 'Mandol' seems to be erroneous (e.g. Bech 1990) 
and therefore Castillejo (1990) assigned eight specimens with an^. gilvus morphology 
from Serra de Pandols near Gandesa (Province of Tarragona) as topotypes [see also 
Castillejo and Rodriguez (1991)]. Subsequently, A. gilvus was redescribed by Garrido 
(1992). Afterwards, the species has also been found in the Provinces of Valencia, Ter- 
uel and Albacete [Borreda (1994), figure 15 in Castillejo (1997), figure 1 in Quinteiro 
et al. (2005)]. Arion gilvus reaches a length of up to 65 mm when extended. It has a 
yellowish to brown dorsum that gets lighter downwards at the sides and dark lateral 
bands that have a yellowish grey line on their upper side (Figure 1). The sole is white or 
evenly yellowish and the mucus is pale yellow (Torres Minguez 1925, Bech 1990, Gar- 
rido 1992, Castillejo 1997). The epiphallus, the pedunculus of the bursa copulatrix, 
and the free oviduct join the atrium on a single line with the pedunculus of the bursa 
copulatrix in the middle, as in A. molinae, but in contrast to the latter, the epiphallus 
is longer than the vas deferens (Torres Minguez 1925, Borreda 1994, Castillejo 1997, 
and figures 26-28 in Garrido et al. 1995). 
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Figure 2. Arion gilvus Torres Minguez, 1925 from Serra de Pandols (Valencia, Spain). A dorsal view 
B lateral view C ventral view. 

As illustrated hj Arion ponsi and^. gilvus, the alleged species-specific genital differ- 
ences among the Iberian species of the A. subfuscus complex are very subtle and little is 
known about their intraspecific variation. Moreover, genital differences among arionid 
taxa do not necessarily imply reproductive isolation (Dreijers et al. 2013). Hence, if 
alleged species-specific phenotypic differences in arionids are to be interpreted under a 
phylogenetic species concept, then their correlation with reproductive isolation should 
be corroborated by molecular data. Molecular markers have been very effective in this 
respect (e.g. Pinceel et al. 2005a, b, Quinteiro et al. 2005, Geenen et al. 2006, Jordaens 
et al. 2010). As such, Quinteiro et al. (2005) investigated the taxonomic affinities of 
Iberian Mesarion species using DNA sequence data. Their analysis of the nuclear ri- 
bosomal internal transcribed spacer 1 region (ITSl) showed a polytomy oi Mesarion 
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species, yet, the analysis of the mitochondrial NADH dehydrogenase I (NDl) gene 
suggested a strongly bootstrap supported group of Iberian Mesarion species with a 
continental-Mediterranean distribution {A. paularensis, A. baeticus, A. urbiae, A. an- 
guloi, A. wiktori, and^. gilvus), and an unsupported group of species with an Atlantic 
distribution {A. lusitanicus, A. nobrei, A. fuligineus, A. hispanicus and A. flagellus). In 
addition, the positions of three Pyrenean species {A. lizarrustii, A. iratii, A. molinae) 
remained unresolved. More specifically, the ND 1 data placed A. gilvus as sister taxon 
of^. urbiae and A anguloi. Quinteiro et al. (2005) did not study individuals from the 
Balearic Islands and thus probably did not include A. ponsi. 

Because DNA sequence data do not only provide phylogenetic information, but can 
also serve as DNA barcodes for species identification (Hebert et al. 2003, 2004), we here 
expand on the work of Quinteiro et al. (2005) by 1) characterizing^, gilvus and^. ponsi 
using mitochondrial COI and 16S rDNA gene fragments, and the larger part of the nu- 
clear ITSl region, 2) exploring the phylogenetic affinities of^. gilvus znAA. ponsi wi^m 
the subgenus Mesarion, and 3) providing diagnostic COI barcodes for both species. 



Material and methods 

Information on the species and specimens included here is provided in Table 1 . In 
total, we screened 45 specimens (Table 1). DNA was extracted from small parts of the 
foot using a NucleoSpin Tissue Kit (Macherey-Nagel, Diiren) following the manufac- 
turer's instructions. PGR reactions were done in 25 |^1 reaction volumes that contained 
1.5 mM MgCl^ in 1 x PGR buffer (Qiagen), 0.2 mM of each dNTR 0.2 \M of each 
primer and 0.5 units of Taq polymerase (Qiagen). A fragment of the mitochondrial 
COI and 16S genes was amplified using primer pairs LCO1490 and HC02198 (Fol- 
mer et al. 1994) and l6Sar and l6Sbr (Palumbi 1996), respectively. The nuclear ITSl 
region (except the ± first 30 bp) was amplified using the primer pair ITSIL and 58C 
(Hillis and Dixon 1991). The PGR profile was an initial denaturation step of 5 min at 
95 °C, followed by 35 cycles of 45 s at 95 °C, 45 s at an annealing temperature of 40 °C 
(COI), 42 °C (16S) or 55 °C (ITSl) and 1.5 min at 72 °C, and ending with a final 
extension step of 5 min at 72 °C. PCR products were purified using the GFX PCR 
DNA Purification Kit (GE Healthcare) following the manufacturer's instructions. Pu- 
rified DNA was diluted in 15 |.il of sterile water. PGR-products were bidirectionally 
sequenced using the ABI PRISM BigDye" Terminator vl. 1 Gycle Sequencing Kit and 
run on a ABI3 130x1 Genetic Analyzer. Sequences were assembled in SeqScape v2.5 
(Life Technologies) and inconsistencies were checked by eye on the chromatogram. Se- 
quences were submitted to GenBank under accession numbers KF305196-KF305225 
for GOI, KF356212-KF356245 for 16S and KF385449-KF385469 for ITSl. These 
datasets were supplemented with DNA sequences from GenBank [including a few spe- 
cies of the other Arion subgenera (Table 1)]. We used those of Carinarion as outgroup. 

Sequences were aligned in GlustalW (Thompson et al. 1994) with default settings 
and without subsequent manual adjustments. In each alignment sequences were trimmed 
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to equal length. The final alignments had a length of 504 bp (COI), 408 bp (16S) and 
587 bp (ITSl), and of 1499 bp after concatenating the three fragments. The COI se- 
quences were translated to amino acid sequences to check for stop codons (but none were 
found). The ITSl sequences were also analysed together with those of Quinteiro et al. 
(2005). In this way we could extend our taxon coverage to A. hispanicus Simroth, 1886, 
A. fuligineus yioK\tt, 1845 and A wo^re/ PoUonera, 1889 (Table 1). Because Quinteiro 
et al. (2005) used other ITSl primers, we had to trim this dataset to a length of 378 bp. 
For each gene fragment, and for the concatenated dataset, we constructed Neighbour- 
Joining (NJ) trees (Saitou and Nei 1987) using the Kimura 2-parameter (K2P) model in 
MEGA v5 (Tamura et al. 201 1) with complete deletion of insertions and deletions (in- 
dels). Branch support was evaluated with 1000 bootstrap replicates (Felsenstein 1985). 
Only bootstrap values > 70% were considered as indicating strong support (Hillis and 
Bull 1993). Uncorrected p-distances (hereafter simply referred to as p-distance) were 
calculated in MEGA v5 (Tamura et al. 2011). For these calculations we considered the 
following Molecular Operational Taxonomic Units (MOTUs): 1) the five 16S rDNA 
clades of A subfuscus (SI to S5) defined by Pinceel et al. (2005a), 2) A anguloi and A 
urhiae ]omt\y as a single MOTU (Backeljau et al. 1994, Quinteiro et al. 2005), 3) A. 
wiktori and A. paularensis jointly as a single MOTU (Backeljau et al. 1 996, Quinteiro et 
al. (2005), and 4) A. lusitanicus from Portugal vs. A. lusitanicus from elsewhere as two dif- 
ferent MOTUs (Davies 1987, Castillejo 1998, Quinteiro et al. 2005). Standard errors of 
mean p-distances among taxa and MOTUs were calculated on 1000 bootstrap replicates. 



Results 

Overall 

The alignments comprized 504 bp for COI (196 variable sites), 408 bp for 16S (121 
sites with alignment gaps, 122 variable sites) and 587 bp for ITSl (277 sites with align- 
ment gaps, 64 variable sites). For the concatenated dataset, there was strong support for 
the subgenera Carinarion, Kobeltia (excluding A. wiktori) znd Arion s.s., and for a clade 
of Arion s.s. + Mesarion (including A. wiktori) (Figure 3). The subgenus Mesarion was 
not monophyletic but consisted of (1) a clade of A. flagellus, A. wiktori, A. paularensis, A. 
baeticus,A. urbiae,A. anguloi, and A. gilvus, (2) two haplotypes of A. lusitanicus (lus-79 
and lus-186) that formed a sister group of Arion s.s. [insofar A. lusitanicus is, of course, 
considered as a member of Mesarion; see e.g. Backeljau (1989)], and (3) a number of 
species/clades among which the relationships were mostly unresolved. Within A. sub- 
fuscus (for which the monophyly was not supported) there were five clades (SI to S5), 
with strong support for (S1,S5),S4) and (S2,S3). The mean p-distance (± SE) among 
the Mesarion OTUs (including A. ponsi and A. gilvus) was 0. 168 ± 0.011 (range: 0. 1 1- 
0.22) for COI, 0.134 ± 0.012 (range: 0.058-0.195) for 16S, and 0.022 ± 0.004 (range: 
0.000-0.048) for ITSl (a minimum distance of zero means that the two sequences only 
difi-ered in a number of indels). The mean p-distances (± SE) excluding A. and A. 
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Figure 3. Neighbour-Joining tree (Kimura 2-parameter model) of a 1499 bp concatenated fragment 
(504 bp of the mitochondrial cytochrome c oxidase subunit I (COI) gene, 408 bp of the mitochondrial 
16S rDNA gene, 587 bp fragment of the nuclear internal transcribed spacer 1 (ITSl) region) for the land 
slug subgenus Mesarion. Bootstrap values > 70% are shown at the nodes. For sample codes see Table 1 . 
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gilvus wcK 0.167 ± 0.011 (range: 0.11-0.22) for COI, 0.130 ± 0.012 (range: 0.058- 
0.195) for 16S, and 0.023 ± 0.004 (range: 0.000-0.048) for ITSl. For the concatenated 
dataset these values were 0.108 ± 0.006 (range: 0.071-0.137) (including A and A 
gilvus) and 0.107 ± 0.006 (range: 0.071-0.137) (excluding^, ponsi and A gilvus). The 
phylogenetic trees inferred from the three gene fragments and from the concatenated 
dataset are shown in Appendix, Supplementary Figures 1-4 and Figure 3, respectively. 



Arion ponsi 

The four individuals oiA. ponsijieXAeA four COI and three 16S haplotypes (Appendix, 
Suppkmentary Figures 1-2), yet two 16S hapfotypes only differed by an indel of two 
base pairs at positions 291-292. For both genes A. molinae showed the smallest p-dis- 
tance with A ponsi (COI: mean p-distance 0.142 ± 0.014; 16S: mean p-distance 0.162 
± 0.019), but a sister species relationship with A. molinae was only well-supported by 
16S. There were three ITSl haplotypes for^. ponsi; one of these had a deletion of a 
poly-T stretch of six base pairs at positions 556-561; the other two differed by a dele- 
tion of a G at position 554. These three ITSl haplotypes of^. ponsi clustered within a 
clade oiA. subfuscus S\-5 , A. lizarrustii, A. molinae, A. iratii and A. transsylvanus {Ap- 
pendix, Supplementary Figure 3). The ITSl analysis with the sequences of Quinteiro 
et al. (2005), placed the single remaining A ponsi haplotype in the same clade (mean 
p-distance with the other taxa of this clade = 0.046 ± 0.004), but without bootstrap 
support (Appendix, Supplementary Figure 4). 

As for 16S, the concatenated tree of the three gene fragments showed a sister spe- 
cies relationship between A. ponsi and A. molinae (Figure 3). 



Arion gilvus 

The three ^. gilvus specimens yielded two COI (one synonymous A-G substitution at 
position 366) and one 16S haplotypes. For both genes the smallest mean p-distances 
were observed relative to^. urbiae2indA. anguloi (COI: mean p-distance = 0.145 ± 
0.013; 16S: mean p-distance = 0.134 ± 0.016). The two ^4. gilvus ITSl haplotypes 
reduced to one when considering the stretch that overlapped with the Quinteiro et 
al. (2005) sequences. In this stretch it differed from that of Quinteiro et al. (2005) 
by a deletion of a T at position 349. Separately, none of the three genes provided 
reliable evidence about the sister group relationships of A. gilvus (Appendix, Sup- 
plementary Figures 1-4). Yet, the concatenated tree showed a well-supported sister 
species relationship between A. gilvus and the A. urbiae I A. anguloi clade (mean 
p-distance = 0.021 ± 0.003) (Figure 3). Mean p-distances within this A. urbiae I 
A. anguloi clade (in which A anguloi was, paraphyletic) were P = 0.041 ± 0.006 for 
COI, P = 0.023 ± 0.006 for 16S, P = 0.004 ± 0.002 for ITSl and P = 0.020 ± 0.003 
for the concatenated dataset. 
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Discussion 

The NJ-tree of the concatenated dataset confirms the major outcomes of previous phy- 
logenetic studies, viz. 1) a strong support for the monophyly of the subgenus Carinarion 
(Geenen et al. 2006), 2) a clade of Arion s.s. and non-Portuguese A lusitanicus (Quin- 
teiro et al. 2005), 3) A. wiktori clustering with Mesarion species, in particular with A. 
paularensis (Quinteiro et al. 2005) instead of with Kobeltia species (Castillejo 1998), and 
4) the strong differentiation within^, subfuscus s.s. that consists of at least, five phyloge- 
netic species (Pinceel et al. 2005a). It therefore seems that the analysis of COI, 16S and 
ITSl DNA sequences yields relevant taxonomic information with respect to the charac- 
terisation of arionid species that have been described under the morphospecies concept. 

Because Arion gilvus and Arion ponsi were originally described on morphological 
grounds they are to be interpreted as morphospecies. This phenetic morphological 
distinction, however, correlates well with a phenetic separation based on mtDNA dis- 
tances. Indeed, the overall mean p-distance among the Mesarion MOTUs (excluding 
A. ponsi and A. gilvus) dealt with in this study is 16.7% for COI and 13% for 16S. 
As such, the mean p-distances between A ponsi and^. molinae (COI: 14.2%; 16S: 
16.2%) or between A. gilvus and A. urbiae (COI: 14.5%; 16S: 13.4%) are perfectly 
comparable with the mean p-distances among the other MOTUs and morphospecies 
in Mesarion. Hence, with respect to mtDNA divergence, both A. ponsi and A. gilvus, 
behave as other Mesarion species or putative species-level MOTUs. 

Obviously, the strong COI differentiation ^mong Mesarion taxa, and of A. ponsi 
and A. gilvus in particular, suggests that DNA barcoding may be a suitable identi- 
fication tool for these animals. Yet, this may be a too simplistic conclusion, since 
stylommatophorans may show extremely high intraspecific mtDNA divergences of 
sometimes up to 27% (K2P-distances, but note the uncorrected p-distances are al- 
most similar) (Thomaz et al. 1996, Chiba 1999). In addition, Davison et al. (2009) 
showed that in the Stylommatophora the mean interspecific K2P-distances (± 3%) 
can be substantially lower than the mean intraspecific K2P-distances (± 12%). Un- 
der these conditions, it becomes very difficult to define generally applicable thresh- 
olds that distinguish between intra- and interspecific sequence divergences. Such 
thresholds are normally associated with DNA barcoding gaps (Hebert et al. 2003), 
but Davison et al. (2009) were unable to detect DNA barcoding gaps in the taxa they 
studied. Nevertheless, Davison et al. (2009) suggested a pragmatic 4% threshold to 
separate intra- and interspecific values, but at the same time they also concluded that 
DNA barcoding in itself is insufficient to identify and/or detect stylommatophoran 
species. Unfortunately, our sample sizes were too small to explore eventual DNA 
barcoding gaps in Mesarion. 

Because DNA barcoding on its own may be unreliable for identifying and detect- 
ing species-level taxa in stylommatophorans, it its necessary to backup this sort of data 
with, amongst others, phylogenetic analyses. As such, our phylogenetic trees of the DNA 
sequence data show that the morphospecies A. ponsi and A gilvus, also represent phylo- 
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genetic species, since both form well-supported clades that are "significantly" associated 
with well-defined, but morphologically different sister species. For A. ponsi, the sister spe- 
cies appears to be A molinae, the distribution range of which is located in NE continental 
Spain (Castillejo 1997), i.e. north of, and facing, the Balearic Islands. Conversely, the 
sister taxon of A gilvus is the "tandem" of A urbiae and A. angubi, two species that have 
been synonymized by Backeljau et al. (1994) and that jointly should be referred to as A. 
urbiae. Our DNA sequence data on COI, 16S and ITSl (e.g. Figure 3), as well as those 
on NDl and ITSl of Quinteiro et al. (2005) are in line with this. As such, the distribu- 
tion range of A urbiae is situated northwest of, and probably adjacent to, that of A. ^Ivus. 
Thus, for both the species pairs A ponsi I A. molinae and A gilvus I A urbiae, the distribu- 
tion ranges appear at least consistent with the suggested sister group relationships. 

In conclusion, the present work shows that A. ponsi and^. gilvus clearly differ from 
A. subfuscus or any other currently recognized arionid species. As such, former records 
of A. subfuscus from Menorca (e.g. GasuU and van Regteren Altena 1970, Mateo 1993, 
Beckmann 2007) almost certainly refer to A. ponsi. Similarly, probably all reports of A. 
subfuscus in the regions of Valencia and Albacete involve A. gilvus (e.g. Borreda 1994, 
Borreda and CoUado 1996). Finally, Borreda (1994) wondered about the eventual re- 
lationship between A. subfuscus from Menorca and A. gilvus. The current data confirm 
unambiguously that these are two different species, with the former being A. ponsi. Yet, 
the overall phylogenetic relationships within Mesarion and many other A. subfuscus- 
like taxa remain to be resolved. In this context, one of the main questions is whether 
Mesarion in its present use is a monophyletic taxon. At the same time one may wonder 
about the relationships with the subgenus Arion s.s., with which Mesarion seems to 
form a well-supported clade (Figure 3). 
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Supplementary Figure I . Neighbour-Joining tree (Kimura 2-parameter model) of a 504 bp fragment 
of the mitochondrial cytochrome c oxidase subunit I (COI) gene for the land slug subgenus Mesarion. 
Bootstrap values > 70% are shown at the nodes. For sample codes see Table 1 . 
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Supplementary Figure 2. Neighbour-Joining tree (Kimura 2-parameter model) of a 408 bp fragment 
of the mitochondrial 16S rDNA gene for the land slug subgenus Mesarion. Bootstrap values > 70% are 
shown at the nodes. For sample codes see Table 1. 
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Supplementary Figure 3. Neighbour-Joining tree (Kimura 2-parameter model) of a 587 bp fragment 
of the nuclear internal transcribed spacer 1 (ITSl) region for the land slug subgenus Mesarion. Bootstrap 
values > 70% are shown at the nodes. For sample codes see Table 1 . 
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Supplementary Figure 4. Neighbour-Joining tree (Kimura 2-parameter model) of a 378 bp fragment 
of the nuclear internal transcribed spacer 1 (ITSl) region for the land slug subgenus Mesarion. This figure 
also includes the Iberian Mesarion ITSl sequences of Quinteiro et al. (2005) Bootstrap values > 70% are 
shown at the nodes. For sample codes see Table 1. 



